## all packages
library(rms)
library(ggpubr)
library(tidyverse)
library(gt)
library(broom)
library(multiwayvcov)
library(dplyr)
library(ggplot2)
library(cowplot)
library(magick)
library(readr)
library(gtsummary)

#### load data ####
data <- readRDS("/data/dataset_leip_dataverse.rds")

#### prep data ####
with_leip <- data %>% 
  filter(abs(x) < 100) %>% 
  mutate(`Z` = pcteligible2013,
         `Z` = `Z` >= median(`Z`))

#### plot ####
f5p1 <- with_leip %>% 
  mutate(`Z` = ifelse(`Z`, "High eligibility", "Low eligibility"),
         `Z` = factor(`Z`, levels = c("Low eligibility", "High eligibility"))) %>% 
  ggplot() +
  geom_point(aes(x, dturnout20142010), alpha = 0.5, color = "grey", size = 0.5) +
  geom_smooth(aes(x, dturnout20142010, group = x < 0), method = "loess", color = "black") +
  facet_wrap(~Z) +
  scale_y_continuous(limits = c(-20, 10)) +
  geom_vline(xintercept = 0) +
  theme_bw() +
  labs(
    title = "Change in Turnout 2014-2010",
    x = "Distance from Nearest State Border",
    y = "Change in Voter Turnout"
  )

f5p2 <- with_leip %>% 
  mutate(`Z` = ifelse(`Z`, "High eligibility", "Low eligibility"),
         `Z` = factor(`Z`, levels = c("Low eligibility", "High eligibility"))) %>% 
  ggplot() +
  geom_point(aes(x, dturnout20162012), alpha = 0.5, color = "grey", size = 0.5) +
  geom_smooth(aes(x, dturnout20162012, group = x < 0), method = "loess", color = "black") +
  facet_wrap(~Z) +
  scale_y_continuous(limits = c(-20, 10)) +
  geom_vline(xintercept = 0) +
  theme_bw() +
  labs(
    title = "Change in Turnout 2016-2012",
    x = "Distance from Nearest State Border",
    y = "Change in Voter Turnout"
  )


## display plot
ggarrange(f5p1, f5p2, nrow = 2)